%this function calculates the ensemble average of Reynolds stresses,
%Swirling Strength, and Vorticity. In addition, it calculates and saves
%instantaneous fields of fluctuations, vorticity, and swirling strength.
%The inputs are:

%pathFiles: Composite array containing the entire path of the vector fields
%folFluc: folder in which the fluctuation fields are going to be stored
%folSwWz: folder in which the instantaneous fields of vorticity and swirling
%   strength are going to be stored
%Uxavrg, Uyavrg, and Uzavrg: ensemble average of velocity components
%windowProcess: Region of interest to process
%FAC: Correction factor for velocity fields

%created by Ricardo Mejia-Alvarez. Urbana, IL. 07/02/09
%based on the function 'ReynoldsStress' (07/30/08)

%modified:  01/03/2010
%           01/25/2010
%           05/30/2010
%           09/19/2010
%           09/23/2010

function [RestressAv,WzAv,WzRMS,SwAv]...
    = ReStressVortSwirl(pathFiles,folFluc,folWz,folSw,Uavrg,windowProcess,FAC,invertCoord)

path = char(pathFiles);
    
for j = 1 : length(path(:,1))

    [~,I,J,Dx,Dy,U] = matrix_Cell( path(j,:) );

    if j == 1 
        SC = size( U );
        X = cat( 3 , cell2mat( U(1) ) * invertCoord(1) , cell2mat( U(2) ) * invertCoord(2) );
    end

    if SC(2) ~= 8
 
        CHC = cell2mat( U(5) );
        U = cat( 3 , cell2mat( U(3) ) * invertCoord(1) ,...
            cell2mat( U(4) ) * invertCoord(2) , zeros( J , I) );

    else
        
        CHC = cell2mat( U(8) );
        U = cat( 3 , cell2mat( U(4) ) * invertCoord(1) ,...
            cell2mat( U(5) ) * invertCoord(2) , cell2mat( U(6) ) * invertCoord(3) );

    end

    %reducing size of matrices
    if j == 1
        [~ , ~ , X] = WindowFileMat(windowProcess , I , J , X);
    end
    [~ , ~ , U] = WindowFileMat(windowProcess , I , J , U);
    [I , J , CHC] = WindowFileMat(windowProcess , I , J , CHC);
    
    if j == 1
        WzAv = curl( X(: , : , 1) , X(: , : , 2) , Uavrg(: , : , 1) , Uavrg(: , : , 2) );
        WzRMS = zeros( J , I);
        SwAv = zeros( J , I);
        RestressAv = zeros(J , I , 6);
        CHCtotal = zeros(J , I);
        CHCtotalRestress = zeros(J , I , 6);
        ti = num2str(I);
        tj = num2str(J);
    end
    
    if FAC ~= 1
        U = U * FAC;
    end
    
    [CHC] = depurationCHC(CHC);   %identifies places with bad vectors as generated by Insight 3G flag
    [rowWrong , colWrong] = find( CHC == 0 );   %finding coordinates for wrong vectors

    CHCu = repmat( CHC , [1 , 1 , 3] );
    CHCRestress = repmat( CHC , [1 , 1 , 6] );

    %substitutes NANs by zeros
    U = GetRidofNaNs(U);
    
    %velocity fluctuations
    ufluc = (U - Uavrg) .* CHCu;    %eliminates wrong vectors as well
    Ufluct = cat(3 , X , ufluc);
    [Ufluct] = mixingMat(Ufluct);
 
    %Reynolds stresses
    Restress = cat(3 , ufluc(: , : , 1).^2 , ufluc(: , : , 2).^2 , ufluc(: , : , 3).^2 ,...
        ufluc(: , : , 1) .* ufluc(: , : , 2) , ufluc(: , : , 1) .* ufluc(: , : , 3) ,...
        ufluc(: , : , 2) .* ufluc(: , : , 3) );
    RestressAv = RestressAv + Restress;

    %replacing zero entries in velocity matrices by median of
    %neighbors. This is used to reduce noise in computation of
    %derivatives.
    Uinstx = substituteMedian(rowWrong , colWrong , U(: , : , 1));
    Uinsty = substituteMedian(rowWrong , colWrong , U(: , : , 2));

    %vorticity of fluctuations
    [VorticityZ] = curl( X(: , : , 1) , X(: , : , 2) , Uinstx , Uinsty );
    VorticityZ = (VorticityZ - WzAv).* CHC;
    WzRMS = WzRMS + VorticityZ.^2;
    Wz = cat(3 , X , VorticityZ);
    [Wz] = mixingMat(Wz);
    Wz = sortrows(Wz,[2,1]);

    %swirling Strength
    [Sw] = SwirlingStrength(Dx , Dy , Uinstx , Uinsty , VorticityZ);
    Sw = Sw .* CHC;
    SwAv = SwAv + Sw;
    Sw = cat(3 , X , Sw);
    [Sw] = mixingMat(Sw);
    Sw = sortrows(Sw,[2,1]);

    %calculating sum of good values
    CHCtotal = CHCtotal + CHC;
    CHCtotalRestress = CHCtotalRestress + CHCRestress;

    %these lines save each individual fluctuations field
    title = strcat('VARIABLES= "X", "Y", "U", "V", "W", ZONE',' I=',ti,' J=',tj,' F=POINT');
    name = filename(path(j,1:length(path(j,:))-4));   %this is the name of the file without extension
    nameP = strcat('\',name,'-fluc.dat'); %this is the new name of the file
    saver(folFluc,nameP,title,Ufluct)

    %these lines save each individual field of vorticity
    %fluctuations
    title = strcat('VARIABLES= "X", "Y", "<greek>W</greek><sub>Z</sub>", ZONE',...
        ' I=',ti,' J=',tj,' F=POINT');
    nameP = strcat('\',name,'-Vort.dat'); %this is the new name of the file
    saver(folWz,nameP,title,Wz)

    %these lines save each individual field of swirling strength
    title = strcat('VARIABLES= "X", "Y", "<greek>l</greek><sub>ci</sub>", ZONE',...
        ' I=',ti,' J=',tj,' F=POINT');
    nameP = strcat('\',name,'-Swirl.dat'); %this is the new name of the file
    saver(folSw,nameP,title,Sw)
            
end
    
RestressAv = RestressAv ./ CHCtotalRestress;
SwAv = SwAv ./ CHCtotal;
WzAv = WzAv ./ CHCtotal;
WzRMS = sqrt( WzRMS ./ CHCtotal );

[CHCtotal] = depurationCHC(CHCtotal);   %identifies places with wrong vectors as generated by Insight 3G flag
[rowWrong , colWrong] = find( CHCtotal == 0 );   %finding coordinates for wrong vectors

%substitutes wrong vectors by zeros
for m = 1 : length(rowWrong)
    RestressAv(rowWrong(m) , colWrong(m) , :) = 0;
    SwAv(rowWrong(m) , colWrong(m) ) = 0;
    WzAv(rowWrong(m) , colWrong(m) ) = 0;
    WzRMS(rowWrong(m) , colWrong(m) ) = 0;
end

%substituting NaNs by zeros
SwAv = GetRidofNaNs(SwAv);
WzAv = GetRidofNaNs(WzAv);
WzRMS = GetRidofNaNs(WzRMS);
RestressAv = GetRidofNaNs(RestressAv);

end         